##### for Extended Data Figure 12 b,c

rm(list=ls())
options(stringsAsFactors = F)
library(Seurat)
library(ggplot2)
setwd("/Users/gin/Documents/onedrive/data/public_database/scRNA/E-GEOD-76312/")

######读取表达矩阵#####---------
sce <- CreateSeuratObject(Read10X('./sce/'),"sce")

######PCA降维#####---------
library(Seurat)

sce <- NormalizeData(object = sce, 
                      normalization.method = "LogNormalize",
                      scale.factor = 10000)

sce <- FindVariableFeatures(sce, selection.method = "vst", nfeatures = 2000) 

sce <- ScaleData(sce,features = NULL)


sce <- RunPCA(sce, features = VariableFeatures(object = sce))
sce <- FindNeighbors(sce, dims = 1:10)
sce <- FindClusters(sce, resolution = 0.1) #resolution越高，细胞群越多，并不是越多越好，需要根据自己的样本特征调整

sce <- RunUMAP(sce, dims = 1:10) 
sce <- RunTSNE(sce, dims = 1:10)


DimPlot(sce, reduction = "umap",label = T)
DimPlot(sce, reduction = "tsne",label = T)
#save(pbmc,file="pbmc_umap.Rdata")



# FGFR1 ENSG00000077782 BMPR2 ENSG00000204217 BMPR1A ENSG00000107779 BMPR1B ENSG00000138696
FeaturePlot(sce,reduction = "tsne",features = c("ENSG00000077782","ENSG00000204217","ENSG00000107779","ENSG00000138696"),order=T)
# FGFR1 ENSG00000077782
# BMPR2 ENSG00000204217
# BMPR1A ENSG00000107779
# BMPR1B ENSG00000138696

### 特定群
ba_sce = sce[,sce@meta.data$seurat_clusters %in% 0]
normal_sce = sce[,sce@meta.data$seurat_clusters %in% 1]



##### 提取表达量 ##### ------
#FGFR1 ENSG00000077782 BMPR2 ENSG00000204217 BMPR1A ENSG00000107779 BMPR1B ENSG00000138696
ba <- as.matrix(ba_sce@assays$RNA@counts)
normal <- as.matrix(normal_sce@assays$RNA@counts)

#meta <- as.matrix(pbmc@meta.data$seurat_clusters)
#meta[meta=="0"]="CML"
#meta[meta=="1"]="Normal"
gene <- c("ENSG00000077782","ENSG00000066468","ENSG00000068078","ENSG00000160867", "ENSG00000107779","ENSG00000138696","ENSG00000204217")
ba=ba[gene,]
rownames(ba)=c("FGFR1","FGFR2","FGFR3","FGFR4","BMPR1A","BMPR1B","BMPR2")
ba <- t(log2(ba+1))
ba = ba[rowSums(ba)>0,]
write.csv(ba,"ba.csv")

normal=normal[gene,]
rownames(normal)=c("FGFR1","FGFR2","FGFR3","FGFR4","BMPR1A","BMPR1B","BMPR2")
normal <- t(log2(normal+1))
normal = normal[rowSums(normal)>0,]
write.csv(normal,"normal.csv")


##### GSEA ##### -----
### data processing
#### deg and id 转换
markers_df <- FindMarkers(object = sce, ident.1 = 0,ident.2 = 1, min.pct = 0.01,logfc.threshold = 0.01,thresh.use=0.99)
library(biomaRt)
ensembl<-useMart("ensembl","hsapiens_gene_ensembl")
attributes<-listAttributes(ensembl)
value<-rownames(markers_df)
attr<-c("hgnc_symbol","ensembl_gene_id")
ids<-getBM(attributes=attr,
           filters="ensembl_gene_id",
           values=value,
           mart=ensembl)
colnames(ids) = c("gene","ensembl")

ids = ids[-13646,]
rownames(ids) = ids$ensembl
deg = markers_df[rownames(ids),]
deg$symbol = ids[,1]
deg = na.omit(deg)

deg = deg[!(deg$symbol==""),]
deg = deg[!duplicated(deg$symbol),]
rownames(deg) = deg$symbol

### GSEA
geneList=deg$avg_log2FC
names(geneList)= rownames(deg)
geneList=sort(geneList,decreasing=T)
head(geneList)
library(ggplot2)
library(clusterProfiler)
library(org.Hs.eg.db)


#选择gmt文件(MigDB中的全部基因集)
gmtfile ='pathway.gmt'
pro1 = "pathway"
pro2='N_vs_CML'
#31120个基因柒
#GSEA分析
library(GSEABase)
geneset <-read.gmt(gmtfile)
length(unique(geneset$term))
egmt <-GSEA(geneList,TERM2GENE=geneset,
            minGSSize=1,pvalueCutoff=0.99,verbose=FALSE)

gsea_results_df <-egmt@result

write.csv(gsea_results_df,file=paste0(pro1,"_",pro2,'_gsea_results_df.csv'))

library(enrichplot)

gseaplot2(egmt,gsea_results_df$ID[1],
          title=gsea_results_df$Description[1],pvalue_table=F)
ggsave(paste0(pro2,'_',gsub('/','-',gsea_results_df$Description[1]),'.pdf'),width = 6,height = 6)

